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ABSTRACT 

Solutions to the energy-independent (gray) radiative transfer equations are compared to results of Monte 
Carlo simulations of the 56 Ni and 56 Co radioactive decay y-ray energy deposition in supernovae. The com- 
parison shows that an effective, purely absorptive, gray opacity, K y ~ (0.06 ± 0.01)7,. cm 2 g _1 , where Y e is the 
total number of electrons per baryon, accurately describes the interaction of y-rays with the cool supernova 
gas and the local y-ray energy deposition within the gas. The nature of the y-ray interaction process 
(dominated by Compton scattering in the relativistic regime) creates a weak dependence of (c y on the optical 
thickness of the (spherically symmetric) supernova atmosphere: The maximum value of tc y applies during opti- 
cally thick conditions when individual y-rays undergo multiple scattering encounters and the lower bound is 
reached at the phase characterized by a total Thomson optical depth to the center of the atmosphere r e <, 1. 
However, the constant asymptotic value, K y = 0.050 Y e cm 2 g _1 , reproduces the thermal light curve due to 
y-ray deposition for Type la supernova models to within 10% for the epoch from maximum light to t = 1200 
days. Our results quantitatively confirm that the quick and efficient solution to the gray transfer problem 
provides an accurate representation of y-ray energy deposition for a broad range of supernova conditions. 
Subject headings: gamma rays: theory — nuclear reactions, nucleosynthesis, abundances — 
radiative transfer — supernovae: general 


1. INTRODUCTION 

Explosive silicon burning at high temperatures and densities 
during supernova explosions can produce unstable isotopes of 
iron group elements (Truran, Arnett, & Cameron 1967; 
Bodansky, Clayton, & Fowler 1968). The decay of these iso- 
topes and subsequent thermalization of the decay products can 
generate the energy needed to power the observed supernova 
optical display (Baade et al. 1956; Colgate & White 1966; 
Colgate & McKee 1969; Weaver & Woosley 1984) especially 
during the slowly evolving later phases. Numerical simulations 
invoking radioactive decay as a source of energy are able to 
reproduce observed supernovae light curves in detail. Not until 
SN 1987A, however, has the radioactive scenario been directly 
confirmed by the detection of y-ray lines from 56 Co decay 
(Matz et al. 1988; Arnett et al. 1989; Palmer et al. 1993) though 
indirect evidence for the fresh synthesis of 56 Ni has been 
observed in optical (Axelrod 1980) and infrared (Varani et al. 
1990) spectra by measurements of the evolution of thermally 
excited Co emission features consistent with 56 Co decay. 

A practical consideration in supernova dynamics modeling 
is determining the local rate of radioactive energy deposition 
(or local heating rate). This rate provides the source term in the 
energy equation as part of the system of equations describing 
the physical state of the atmosphere. It can also be used 
directly in certain approximations for estimating rates of ion- 
ization and excitation by energetic electrons in statistical equi- 
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librium (Axelrod 1980; Chugai 1987; Swartz 1991; Xu & 
McCray 1991; Kozma & Fransson 1992) and nonequilibrium 
(Fransson & Kozma 1993) equations. The radioactive source 
distribution can be provided by nucleosynthetic yields from 
explosion simulations, but the local heating rate is determined 
by the rate at which y-rays and positrons deposit energy at 
various points in the gas as they travel through the atmo- 
sphere. (It is customary to assume that positrons will annihi- 
late in the immediate vicinity of where they are produced. This 
is justified in terms of the relatively shorter range of electrons 
and positrons when compared with y-rays and the likelihood 
of a tangled magnetic field sufficiently strong to inhibit the 
diffusion of charged particles, though see below.) Thus, in prin- 
ciple, a set of y-ray transport equations must be solved simulta- 
neously with the equations of hydrodynamics and equations of 
state in order to estimate observable quantities from models of 
supernovae. (The hydrodynamical evolution is often simplified 
by assuming homologous expansion which is justifiable once 
the ejecta have increased in radius by a factor of ~ 10 over the 
presupernova stellar size.) 

Examples of the various methods for treating the y-ray 
transport and energy deposition include Monte Carlo tech- 
niques (Colgate, Petschek, & Kriese 1980; Ambwani & Suther- 
land 1988, hereafter AS; The, Burrows, & Bussard 1990), 
multigroup (Weaver, Axelrod, & Woosley 1980), and gray 
(Sutherland & Wheeler 1984, hereafter SW) radiative transfer 
methods as well as analytic approximations (Colgate et al. 
1980; Weaver et al. 1980; Arnett 1979; Ensman & Woosley 
1988). It has been argued (SW) that the complex, multiple- 
scattering character of y-ray radiative transfer in supernova 
envelopes (where Compton scattering dominates) can, for the 
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purposes of calculating energy deposition, be approximated by 
purely absorptive transfer. The plausibility of this assumption 
hinges on the following features of Compton scattering: (i) the 
scattering cross section is forward-peaked, but forward- 
scattering leads to little energy transfer to the gas, and (ii) when 
a large-angle scattering takes place, there is significant energy 
transfer: 

E' y - E y /[1 + (£ y /m c c 2 )( 1 — cos 0)] , 

where E y is the initial photon energy, E' y is the energy after the 
photon scatters from an electron (mass m e ) at rest through 
angle 9. Thus a ~2 MeV y-ray (a representative energy for 
56 Ni and 56 Co decay) gives up >0.8 of its energy if it scatters 
through 6 > 90°. Crudely, then, one can think of the photon as 
proceeding along a nearly rectilinear path until it suffers a 
large-angle scattering, thereby giving up most of its energy. 
Ultimately, the credibility of the pure absorption approach 
depends on its producing results that are reasonably close to 
those from more realistic (i.e., Monte Carlo) simulations. 

The purpose of the present work is to compare the results 
from the physically accurate but computationally intensive 
Monte Carlo simulations to solutions of transfer equations 
assuming a purely absorptive, energy-independent, opacity. 
Section 2 outlines the Monte Carlo and radiative transfer 
methods and the best fit value of the y-ray opacity is deter- 
mined in § 3. We find the purely absorptive effective opacity, 
K y , can be expressed as K y = Y e , where Y e is the (total) 
number of electrons per baryon. In this form, K y is independent 
of the radioactive source distribution, mass density, and ele- 
mental abundance distributions (except for the dependence on 
abundances through Y e ). We find a weak dependence of K y on 
the total electron scattering optical depth, x e , that results from 
a transition from a multiple to a single scattering environment 
encountered by the y-rays in the Monte Carlo simulations. The 
best fit value of K y decreases from K y ~ (0.065 + 0.005) Y e 
cm 2 g - 1 at x e > 10 to K y ~ 0.050 Y e cm 2 g“ 1 asymptotically for 
x e <, 1. A constant value for K y of 0.050 F ( , reproduces the 
Monte Carlo deposition and Type la light curves to within 
10%. 

By combining the effective opacity with a known radioactive 
source distribution in a supernova dynamic model, the rele- 
vant transfer equation becomes a linear first-order differential 
equation with constant coefficients and can be reduced to 
quadrature. Similarly, the Schwarzschild-Milne equations can 
be used to obtain the moments of the radiation field and hence 
the local deposition rate. In most practical applications, 
however, numerical integrations are required because the 
source function cannot be readily expressed analytically. Our 
radiative transfer algorithm to determine the local deposition 
function from known mass density and source distributions 
in spherical symmetry is available upon request to 
pgs@snowdrop.physics.mcmaster.ca. The processing time 
requirement for this algorithm (typically ~0.05 CPU seconds 
for a workstation for the full deposition function calculation) is 
a small fraction (~10~ 5 ) of the time required for a typical 
Monte Carlo simulation following 500,000 decay events. 

2. NUMERICAL METHODS 

The most abundant radioactive isotope produced in super- 
novae is 56 Ni which decays with a 6.1 day half-life to 56 Co 
which in turn decays with a 77.7 day half-life to 56 Fe. The 
decays produce energetic y-rays and positrons with initial ener- 
gies 0.158 to 3.640 MeV (Lederer & Shirley 1978; AS). Super- 


novae yield other, less abundant, radioactive isotopes with 
longer life-times including 57 Co, 44 Ti, and 22 Na. Typical 
photon energies for these sources range from < 100 keV to 
~1.3 MeV. These long-lived isotopes become important only 
at very late times, t > 1000 days after the explosion (e.g., 
Woosley, Pinto, & Hartmann 1989), and are not considered in 
this work. One of the more unique aspects of the radioactive 
decay model for supernovae is that the ambient plasma 
remains cool (T e ~ 1 eV) and nearly neutral while the primary 
heat source is characterized by photons that are six orders of 
magnitude more energetic. This can be seen by equating the 
energy density due to the slow radioactive energy release to 
that of an equivalent blackbody: Assuming the gas is optically 
thick to y-rays then the mean intensity, J ~ e rad /4ruc y , as shown 
below, where e rad ~ 7 x 10 9 exp ( — t/r Co ) ergs s -1 g~ 1 is the 
generation rate per gram at times t P r Ni (r Ni = 7.58 x 10 s s 
and r Co = 9.80 x 10 6 s are the 56 Ni and 56 Co decay times 
[AS]) and K y ~ 0.025 cm 2 g“ This produces an energy 
density e = 4nJ/c ~ e rad /cic y equivalent to that of a blackbody 
at T ~ (itJ/a) 114 ~ 6000 exp ( — t/Qxc) K, where a is the 
Stefan-Boltzmann constant. 

The decay photons primarily interact with the surrounding 
gas through Compton scattering. A fraction of the photon 
energy is lost per scattering and the down scattered photon is 
typically destroyed in X-shell photoionization in less than five 
scatterings (AS). These interactions are independent of the tem- 
perature and ionization state of the gas at the low-temperature 
near-neutral conditions typical of supernovae. The decay posi- 
trons and recoil electrons are the primary ionizing and heating 
agents in the process. Their range is only a fraction of the 
distance across the atmosphere and their energy loss is typi- 
cally treated in situ while the y-rays can travel long distances 
and often escape the atmosphere altogether. Thus the problem 
of calculating the local energy deposition rate is equivalent to 
calculating the rate of production of electrons and positrons 
and their initial spectrum. Positron decays account for 19% of 
all 56 Co decays. The average decay positron energy is ~0.66 
MeV (0.125 MeV per decay or ~4% of the total decay energy) 
and the maximum decay energy is 2.46 MeV. This kinetic 
energy is not accounted for in the Monte Carlo energy deposi- 
tion calculations although the two annihilation photons, each 
of energy E y = m e c 2 , are followed. If the kinetic energy is 
deposited very near the site of the decay, then a small addi- 
tional energy deposition term, proportional to the local mass 
fraction of radioactive 56 Ni, should be added to the results 
from either the Monte Carlo calculation or the gray absorp- 
tion calculations which follow. If the assumption that charged 
particles are trapped locally is invalid, as has been suggested 
again recently by Chan & Lingenfelter (1993), then neither the 
Monte Carlo nor gray absorption calculations considered here 
are appropriate. Chan & Lingenfelter (1993) argue that 
radially-combed magnetic fields may facilitate the escape of the 
positrons and other charged particles. But then at most epochs 
a very large fraction of the radioactive decay energy will escape 
and this would appear to be in conflict with our understanding 
of any supernovae whose light curves we believe to be powered 
by radioactivity. 

2.1. Monte Carlo Simulations 

The Monte Carlo code used in this work is essentially that 
described by AS and by Sutherland (1990). This code was orig- 
inally developed to calculate the emergent y-ray and X-ray 
spectra and incorporated Doppler effects associated with the 
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propagation of the photons through a homologously expand- 
ing supernova atmosphere. The present version includes the 
complete 56 Ni and 56 Co decay source spectrum (AS) and inter- 
actions with the gas through Compton scattering, pair pro- 
duction, and photoelectric absorption. Compton scattering is 
described by the Klein-Nishina cross-section and all electrons, 
bound and ionized, are included. Cross sections for the other 
processes are given in AS with the following correction: The 
pair production cross section, <r pair , per electron is 



where 


{ 0 E y < 2 m e c 2 

0.10063(£ y - 1.022) E y < 1.5 MeV (2) 
0.0481 + 0.301(£ y - 1.5) E y > 1.5 MeV 

in units of 10" 27 cm 2 per atom. Here, E y is the photon energy 
in MeV, n e is the total electron density, and n i is the number 
density of element i with nuclear charge Z,. (Compare eq. [2] 
of AS.) 

The prescription for photoelectric absorption in this code is 
the following: 



per electron where <r p °° is the photoelectric cross-section (per 
atom) for the ith element evaluated at 100 keV and the sum is 
over all ion species. The cross-section data are from the tables 
of Veigele (1973). Numerical values for <r pe for the several com- 
positions used in this work are given in Table 1. These cross 
sections are for neutral atoms but are dominated by the effects 

TABLE i 


Model Composition by Mass Fraction* 


Element 

W7 b 

10H' 

Solar 4 

‘H 


3.1 (-1) 

7.7 (-1) 

2 He 

8.2 (-3) 

5.1 (-1) 

2.1 (-1) 

6 C 

9.3 (-5) 

1.5 (-2) 

3.8 (-3) 

7 N 


2.5 (-4) 

9.3 (-4) 

8 o 

1.0 (-1) 

9.1 (-2) 

8.5 (-3) 

l0 Ne 

3.1 (-4) 

3.9 (-3) 

1.5 (-3) 

12 Mg 

1.9 (-2) 

1-1 (—3) 

7.4 (-4) 

I3 A1 

6.0 (-4) 


6.6 (-5) 

,4 Si 

1.3 (-1) 

1.8 (-2) 

8.1 (-4) 

,5 P 

1.5 (-4) 


5.8 (-6) 

16 S 

6.6 (-2) 

1.1 (-2) 

4.6 (-4) 

17 C1 

1.3 (-4) 


4.8 (-6) 

l8 Ar 

1.4 (-2) 

2.0 (-3) 

1.2 (-4) 

* 9 K 

8.1 (-5) 


3.9 (-6) 

20 Ca 

1.1 (-2) 

1.2 (-3) 

7.2 (-5) 

22 Ti 

2.7 (-5) 

2.5 (-5) 

3.3 (-6) 

24 Cr 

7.0 (-4) 

3.4 (-5) 

1.9 (-5) 

25 Mn 

1.3 (-4) 


1.5 (-5) 

26 Fe 

8.0 (-2) 

5.2 (-4) 

1.5 (-3) 

27 Co 

4.1 (-3) 


3.7 (-6) 

28 Ni 

5.7 (-1) 

3.1 (-2) 

8.1 (-5) 

n 

0.50 

0.66 

0.89 

< 

6.1 (-1)' 

2.6 (-2) 

8.1 (-4) 


* Power oflO exponent in parentheses. 
b From Thielemann et al. 1986 model W7. 
' From Woosley 1988 model 10H. 
d From Cameron 1982. 

'In units of 10 “ 24 cm 2 per electron. 


of the innermost electrons which are not expected to be 
ionized, except for H, in the conditions of interest for super- 
novae ejecta. The scaling £“ 3 for photoelectric absorption is 
an excellent approximation (especially between 10 and 1000 
keV), but in any case the details are less important for deposi- 
tion calculations because by the time a y-ray is likely to suffer 
photoelectric absorption its energy has been so degraded (to 
<30 keV) by repeated Compton scattering that the impact 
upon energy deposition is small. (Gn the other hand, in calcu- 
lating emergent X-ray spectra, it is important to describe accu- 
rately the processes of photoelectric absorption and 
fluorescence.) 

The energy deposited locally by any of the above processes is 
then: (i) the energy difference between the incident and scat- 
tered photon, or (ii) the energy difference between the incident 
photon and 2 m e c 2 if a pair is produced (with the creation of 
two 511 keV y-rays), or (iii) the full energy of the photon if it is 
photoelectrically absorbed. The local deposition function is 
defined as the ratio of the local energy deposition rate per unit 
mass to the rate of energy generation per unit mass of radioac- 
tive material. For a simulation following the fate of N decay 
events, the value of the deposition function in the fcth 
(Lagrangian) mass shell is the product of (i) the ratio of the 
total energy deposited in the shell to the total energy of decay 
photons (£ y Nf y E y , where photons of energy E y are distributed 
according to the probability of emission,/,,, as tabulated in, 
e.g., AS) and (ii) the weight M rad /AM t , where M rad is the total 
mass of initially radioactive material and AM* is the mass of 
the shell. This definition has the property that at early epochs, 
when the ejecta are optically thick and essentially all the y-rays 
are trapped near where they are emitted, the local deposition 
function equals the local mass fraction of radioactive material. 

Typical Monte Carlo simulations involve a minimum of 
500,000 decay events and models with 64 zones. Comparison 
with simulations with even more decays confirmed that this 
number of decays ensures errors less than 5% for all zones for 
all models. 

A separate code that assumes a purely absorptive interaction 
of photons with matter was also implemented, and in this case 
the atmosphere was treated as static so that comparison could 
be directly made with the gray transfer approach of the next 
section. 

2.2. Gray Radiative T ransfer 

The transfer equation in spherical symmetry is cast in a 
difference form and integrated along impact parameters paral- 
lel to an observer’s line of sight. We assume the y- ray opacity, 
K y , is independent of energy, £, and purely absorptive so that 
the transfer equation for the energy-integrated intensity, / = 
Jo / £ dE, along a ray is 

5 /* 

±~^- = n-KyPl (4) 

for the incoming (/") and outgoing (/ + ) directions, where z 
denotes the position along the ray, p is the mass density, and 
t] = $oi 1 E dE is the local total y-ray emissivity: >/ = 
/rad e rad p/^ n - Her e, / rad is the initial mass fraction of 56 Ni and 
e rad denotes the time-dependent rate of energy release per gram 
of radioactive material. In cgs units, 

e rad = 3.9 x 10 10 exp (— t/r Ni ) + 6.78 x 10 9 

x [exp ( - f/i Co ) - exp ( - t/t Ni )] . (5) 
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By introducing the optical depth along the ray, 
dx = —K y pdz, 

and defining 


/' = (4ttk : y /e ni )I , 
equation (4) can be written 


dr 

dx 


V “/rad 


( 6 ) 

(7) 

( 8 ) 


Since no radiation is incident from outside the ejecta, I~ = 0 at 
the upper boundary and, from symmetry at z = 0, the lower 
boundary condition is I + = l~ . Integration of equation (8) 
yields 


I'( t.) = nr i+1 )e-^ + 



(9) 


where x t and T i+1 represent optical depth points along the ray 
and At = (t i+1 — t,). Since / rad is constant within each shell, 
equation (9) can be integrated exactly. The impact parameter 
grid is constructed so that all rays are tangent to radial shells 
with a single “ core ray ” passing through the center. 

The radiative transfer code has been tested successfully 
against another transfer code which casts the transfer equation 
in a second-order difference form and uses a Feautrier solution 
along impact parameters. From this solution, variable Edding- 
ton factors can be evaluated and used in a combined moment 
equation constructed using Auer’s transformation, again 
employing a Feautrier solution. This latter code has been 
extensively tested to reproduce various analytic results includ- 
ing the proper asymptotic results in the gray case (e.g., 
Hummer & Rybicki 1971). 

The rate at which energy is deposited locally (in ergs cm “ 3 * * * * * 
s ” *) is 47 iK y pj, where J is the mean intensity : 


trie shells of radii r k , k = 1, ND, the quality of fit, Q, is defined 
to be 


ND 


Q= I 

k= 1 



(12) 


where d and d e are the local values of the deposition function 
calculated using the Monte Carlo and the gray radiative trans- 
fer methods, respectively. The uncertainties a k are those calcu- 
lated per zone in the Monte Carlo calculation and in effect 
reflect “ counting statistics ” (they are obtained by accumulat- 
ing the variance for each zone in the energy deposited by those 
interactions that occurred in each zone — a k is related to the 
“error in the mean” for the deposited energy). We include a 
weight factor w k for each zone to compensate for the following 
(if w k = 1 then Q is essentially a Xnd statistic): First, we wish to 
downplay low mass zones that were used in the zoning to 
achieve resolution for either the density or velocity profiles. 
Secondly, at times when it is of most interest to calculate the 
deposition function, the thermal diffusion timescale in the 
ejecta is short compared with the dynamic timescale. Then the 
luminosity is given by the instantaneous balance between y-ray 
heating and radiative cooling: L = Y. k AM k d k e rad so that 
regions of small d do not contribute significantly to the 
observed luminosity. Thus Q should be weighted to emphasize 
those zones which contribute most to the luminosity. Thus we 
have chosen the weights w k = The gray deposition, 

d a k , is a function of the single parameter, »c y . Therefore, the 
optimum value of K y is determined for each Monte Carlo simu- 
lation by minimizing Q in equation (12). 

For purposes of exposition reference will be made to the net 
deposition (from which the luminosity, L, follows trivially) 
defined as 


D = v 4c AM * 

” M rad 


(13) 


J = -^j>ldm, (10) 

which is the specific intensity integrated over all solid angles 
dco. Thus the deposition function is: 

d _ 4nK y pJ _ 4nK y J (u) 

^ra(l P ^rad 

The deposition function is (compare eqs. [11] and [7]) then 
simply J', the mean of 

The deposition function can be calculated using the spher- 
ically symmetric gray radiative transfer equations once the 
mass density p, initial 56 Ni mass fraction / rad , and the y-ray 
opacity K y are specified on a radial grid. The calculation does 
not require knowledge of the source spectrum but only of the 
total rate of emission, e rad . 

3. RESULTS 

The Monte Carlo solution is the standard by which the 

deposition function calculated using the gray radiative transfer 

method is to be compared. The deposition function is deter- 

mined by the density profile and the distribution of the radio- 

active material and is therefore position dependent. A global 
measure describing the goodness of fit of the gray results to the 

Monte Carlo results must take into account this spatial depen- 
dence. Assuming the atmosphere is composed of ND concen- 


where M rad = l, k AM t / rad t . The net deposition has the pro- 
perty that D y -*■ 1 as x e -* oo and D y ~* 0 as r e -► 0. 

Simulations were performed for a range of plausible super- 
nova conditions including various mass and source distribu- 
tions, (homogeneous) elemental abundances, and y-ray optical 
depths. Two basic models and their time evolution were con- 
sidered. The first model is based on the incinerated white dwarf 
model W7 of Nomoto, Thielemann, & Yokoi (1984) that is 
known to provide a good fit to the observed light curves 
(Nomoto 1984) and spectra (Harkness 1991) of typical Type la 
supernovae. The total ejected mass for this model is ~ 1.4 M s , 
there is ~0.6 M 0 of radioactive material within the inner ~ 1 
M 0 with a peak in this 56 Ni distribution near M r ~ 0.7 (where 
M r is the mass interior to the point r). The central regions of 
the model, M r < 0.1 M 0 , are depleted of radioactive material 
during neutron-rich freezeout, and the outer ~0.4 M 0 is com- 
posed of partially incinerated, intermediate mass elements that 
result from carbon-burning. The detailed abundances are 
replaced with mass-weighted mean values in the present work, 
but the original 56 Ni distribution is retained intact. 

The second model assumes a power-law density profile 
(p oc r -1 *) with uniform abundances. We examined models with 
total masses in the range 1 to 5 M 0 and power law exponents, 
a, in the range 0 (uniform sphere) to 5 with outer radius R = 
v ma Jt, v ma% ~ 1.0 x 10 9 cm s _l . Both a uniform 56 Ni distribu- 
tion and a central source distribution were considered in the 
power law models. For the central source, the radioactive 
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source distribution is represented by a step function with the 
inner 1% of the total mass assumed radioactive. Chemical 
compositions for model W7 and for two representative mix- 
tures are listed in Table 1. The solar composition represents 
conditions found in the outer layers of typical Type II super- 
novae and the mixture 10H represents Type II supernovae 
with all ejected material (core and envelope) homogenized. 
Each model is initialized using ND = 64 Lagrangian mass 
shells and scaled homologously in time. 

3.1. The Effective y-ray Opacity, K y 

It is instructive to determine the effective opacity for mono- 
energetic sources before considering the full 56 Ni and 56 Co 
decay spectrum. The cross sections for the interaction pro- 
cesses included in the Monte Carlo simulations are dominated 
by Compton scattering for photon energies spanning the 56 Ni 
and “Co decay spectral range, 0.158 < E y < 3.640 MeV (Fig. 
1). The (angle-averaged) Klein-Nishina cross section changes 
monotonically by a factor of 4 within this energy band. The 
effective opacity encountered by a decay photon should there- 
fore depend upon its initial energy. Figure 1 shows the best fit 
to the effective opacity (obtained by minimizing Q) for mono- 
energetic 7 -rays as a function of the initial photon energy for 
model W7 and for several power-law models at an interme- 
diate stage of evolution, D y ~ 0.3 and r e > 3 (see below). The 



Fig. 1. — Best-fit effective opacity to mono-energetic sources for several 
models is shown against the initial photon energy. The upper solid line rep- 
resents the total opacity which is composed of photoabsorption (dotted line), 
Compton scattering ( short-dashed line), and pair production ( long-dashed line). 
All models assume the W7 abundances (T e = 0.5) of Table 1, a denotes the 
power-law exponent for power law models and model W7 is the white dwarf 
model of Nomoto et al. (1984). The total x t to the center of the ejecta for these 
models ranges from ~3 to 12. The lower portion of the figure schematically 
illustrates the 56 Ni and 56 Co decay spectra. The height of the vertical lines 
depicts the relative probabilities of the individual decay paths hence the 
strongest lines are the 56 Ni 158 keV(/ 7 = 1.0) and 56 Co 847 keV(/ 7 = 0.9998) 
lines. 


composition for all models is the mass-averaged W7 model 
abundances. The effective opacity is seen to be independent of 
the model density profile and source distribution but varies 
approximately in proportion to the total cross section at E y for 
initial photon energies between ~0.2 and ~6 MeV. The effec- 
tive opacity is between 0.5 and 0.7 of the total opacity in this 
energy range. The effect of the purely absorptive photoioniza- 
tion process causes this fraction to rise at lower energies with 
the effective opacity approaching the photoionization opacity 
at energies E y <j 0.1 MeV, depending on the relative contribu- 
tion of photoionization to the total opacity. 

The Monte Carlo deposition function is shown for several 
initial y-ray energies against the Lagrangian mass, Af r , in 
Figure 2 for a central source in a 1 M Q uniform density model. 
This illustrates the increasing opacity for lower initial energy 
y-rays from another perspective. The deposition functions for 
initial energies E y > 2 MeV are nearly identical reflecting the 
similar effective opacities in this energy range. At lower initial 
photon energies, the deposition tends to accumulate near the 
point of emission reflecting a relatively higher opacity. Since 
each interaction either destroys the photon or emits a scattered 
photon at a reduced energy, the opacity increases following 
every scattering event. Thus the effective opacity, which reflects 
the changing opacity seen by the photon as it scatters and loses 
energy, is larger than the absorptive opacity experienced only 
by photons at a single initial energy. Further, a photon that 
has lost a significant amount of energy in a scattering event will 
have a high probability of another scattering or absorption 
event nearby. 

The fact that the effective opacity depends on the initial 
photon energy suggests that the effective opacity will be time 
dependent for two reasons. First, the photon energy distribu- 
tion changes from a “Ni-dominated spectrum to a “Co- 
dominated spectrum as the supernova evolves (Fig. 1 and 



Fig. 2. — Logarithm of the deposition function for mono-energetic sources 
is shown for a point source (/ r , d = 1 for M, < 0.1 M 0 ) in a constant density, 1 
M 0 , model at t = 30 days (t e = 15). From top to bottom at Af r = 0, the 
deposition corresponds to initial y-ray energies £ r = 0.15, 0.25, 0.5, 1.0, 3.0, 5.0, 
and 10.0 MeV. The effective opacity is highest for the lowest initial photon 
energies, and the deposition is most strongly concentrated near the point of 
origin as a consequence. The effective opacity is nearly constant for 2 
MeV. This is reflected by nearly identical deposition functions for this energy 
range. The W7 model abundances ( Y f — 0.5) were used. 
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Table 1 of AS). The effective opacity for a pure 56 Ni source can 
be expected to differ from that for a pure 56 Co source (and 
from that for a 57 Co or other radioactive source). These differ- 
ences are, however, negligible compared to the effect due to 
expansion of the atmosphere which reduces the optical thick- 
ness and hence the probability that a photon will interact with 
the gas. 

The time dependence (represented by the total electron scat- 
tering optical depth) of the best fit effective opacity for a mono- 
energetic source is shown in Figure 3 for three of the models 
from Figure 1 (Y^ = 0.5). The effective opacity evolves with 
changing optical depth but is seen to be independent of the 
model source and mass distributions. (The error bars indicate 1 
a confidence intervals for k v Since these errors are determined 
by the number of decay events followed in the Monte Carlo 
simulations, they are shown only to indicate the relative errors 
and are chosen conservatively for clarity.) A general trend from 
large K y at high x e to small K y at low x e is evident. Above x e ~ 1, 
K y has a weak power law dependence on x e . At the highest 
optical depths the destruction length is short and nearly all the 
y-rays are thermalized near their points of origin. In this limit, 
x e -* oo and D y -* 1, the deposition functions are nearly identi- 
cal for a range of K y values above ~0.03. This is also evident 
from equations (9)— (1 1), where, in the optically thick limit, 
d -» / rad independent of the value of K r In the opposite limit, 
x e < 1, K y -* 0.028 cm 2 g ^ 1 for photons with initial energy E y = 
1 MeV. This dependence on x e can be understood qualitatively 
as follows: Suppose at some epoch for which ~ 1 (ij is the 
center to surface optical depth for the scattering of, for 
example, a 1 MeV photon through 90 degrees or more) the best 
fit between the Monte Carlo and gray calculations is achieved 
for K t — Kj. At this epoch, some gamma-rays escape without 
depositing energy, and in the Monte Carlo calculation even 
those photons that scatter will not give up all their energy. The 
combination of these two effects will determine k°. However, at 



Fig. 3. — Increase in the effective opacity with increasing total optical depth 
is shown for several models. The models are homologously expanded giving 
rise to a total Thomson optical depth to the center of the spherically symmetric 
atmosphere which scales as x e oc ( “ 2 . The effective opacity is shown for a 1 
MeV source. As x e -* 0, k -» 0.028 cm 2 g” 1 (Y e = 0.5). Error bars represent the 
1 a uncertainty in the best-fit effective opacity. As r e -» oo, the deposition 
function (and hence Q defined by eq.[12]) is nearly identical for a range of 
opacities leading to formally large uncertainties at high optical depth. 


earlier epochs, when tj > 1 very few photons in the Monte 
Carlo calculation will escape, and by virtue of repeated scat- 
terings will give up virtually all their energy. Thus Kj will have 
to be increased to compensate. Conversely, at epochs for which 
x i <§ 1, the “real” Monte Carlo photons rarely scatter more 
than once, almost never give up their entire energy, and K t 
must accordingly be reduced below k° Put more succinctly, 
the effect of multiple scatterings is to increase the “ efficiency ” 
of energy deposition and therefore at earlier epochs k, must be 
slightly increased. That k, approaches 0.028 cm 2 g _ 1 as x e -» 0 
can be confirmed for a central point source with the following 
simple numerical calculation: The effective cross section for 
single-scattering energy deposition by a y-ray line of energy 
0.2 < E y < 10 MeV is dominated by the energy-loss-weighted 
Klein-Nishina cross-section : 


a tt((Ey) 




^^A EdQ 
dtl 


(14) 


where the energy loss is given by 


AE/E y = (E y /m e c 2 )( 1 - cos 0)/[ 1 + (EJm e c 2 )(\ - cos 0)] . 


For E y = 1 MeV, a eff = 0.14 in units of the Thomson cross- 
section (K eff = 0.028 cm 2 g~ 1 for Y e = 0.5). This result has been 
confirmed for an optically thin sphere with central point 
source, where only single-scattering is important and the 
optical depth to the surface is isotropic. This argument must 
also extend to all optically thin situations. Figure 4 displays 
ff eff for a range of y-ray energies based on equation (14). Figure 
4 also shows that the angle-averaged fractional energy loss is 



Fig. 4. — Effective cross section, <r eff , for single-scattering energy deposition, 
equation 14, is shown along with the angle-averaged Klein-Nishina cross 
section, <r KN , and the angle-averaged fractional energy loss, AE/Ey for a range 
of incident y-ray energies. (u Th is the Thomson cross section. The opacity is 
simply related to the Compton scattering cross sections as k = crYJm.) 
Symbols denote the effective opacity for the oc = 0 uniformly distributed source 
model from Fig. 1. For this model, the effective opacity is slightly higher than 
the analytic solution at all energies due to multiple scattering (r f ~ 4 for this 
model). The larger effective opacity at lower energies is the result of the added 
opacity due to photoionization. 
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~0.5 for E y > 1.0 MeV indicating that about half the photon 
energy is lost per scattering, on average, for photons in this 
energy range. Thus, if only one scattering is encountered by a 
photon, we can expect an effective opacity K y ~ 0.5k kn = 
0.5(<t kn YJm), in agreement with Figure 1. This does not 
account for photoionization which adds considerably to the 
effective opacity at low energies. The true effective opacity, K y , 
which accounts for all interaction processes and all subsequent 
scatterings, exceeds the “ single-scattering ” effective Compton 
opacity, k c(( . Thus, k c{[ is an approximate lower limit for K y in 
that K y -* K err in the optically thin limit and for photon ener- 
gies, E y , far above the photoionization threshold. This reason- 
ing also implies, as confirmed by the Monte Carlo simulations, 
that the largest energy loss occurs during the first scattering 
event. 

These arguments extend to the full 56 Ni and 56 Co decay 
spectrum. The time evolution of the effective opacity for several 
models is shown in Figure 5 using the full time-dependent 
decay spectrum. The trends noted above are reproduced by the 
full spectrum: the best-fit values of K y are independent of the 
source and density distributions and K y has a weak dependence 
on T e , decreasing from ~ 0.035 ± 0.01 at r e > 3 to ~ 0.025 for 
t„ < 1. The effective pure absorption cross section for single 
scattering (eq. [14]) is now: 


ff abs — 


ZyfyEyO'ffjEy) 
■t-yfy Ey 


(15) 


where the sum is taken over the line spectrum (f y is the prob- 
ability for a given line). The result of doing the above calcu- 
lation for the line spectrum of 56 Co decay is <r abs ~ 0.128 in 
units of the Thomson cross section (x y = 0.05 1 Y e cm 2 g“ *). 

Figure 6 shows that K y for different compositions can be 
written as K y = K y Y e , where Y e is the number of electrons per 
baryon (in all the preceding figures, Y e = 0.5). This is the well- 
known Compton scattering dependence on the number of elec- 
trons: For Compton scattering, the effective opacity can be 
written as K y = (n e o KN /p) where the total (bound and free) elec- 
tron density is n e = (p/m)Y e and m denotes the atomic mass 
unit. Thus K y can be expressed as tc y = ( <r/m)Y e or, equivalently. 



Fig. 5. — Evolving effective opacity for several models is shown as in Fig. 3. 
The full 56 Ni and 56 Co decay spectrum is used here instead of a mono- 
energetic source as in Fig. 3. 



Fig. 6. — Evolving effective opacity for models with different abundances 
(Table 1) is shown as in Fig. 3. In this figure, k 0 = kJY, is plotted while in 
previous figures k t was plotted with Y t = 0.5. The effective opacity, k v depends 
on the abundances, to a very good approximation, only through the Compton 
scattering dependence on the number of available electrons per unit mass. 


as K y = K y Y e . The values of K y shown in Figure 6 and in pre- 
vious figures are consistent with values quoted in the literature. 
For example, Colgate et al. (1980) found K y = 0.028 cm 2 g _l 
for a pure Ni 0.5 M e uniform sphere Type la supernovae 
model. Woosley et al. (1989) quote a value of 0.05y c for the 10H 
composition model for SN 1987A at late times. 

3.2. Deposition Functions and Light Curves 

The quantity Q provides a global statistical measure of the 
agreement of the gray deposition function to the results of the 
Monte Carlo simulation. By minimizing Q, the best-fit value of 
the parameter x y is determined. What is also needed is a 
measure of the accuracy of the results parameterized by x r 
Two physical properties that characterize the y-ray transport 
phenomenon and that are fundamentally relevant to super- 
nova studies are the local heating rate and the light curve. 
These are related to the local deposition function, d k , and the 
net deposition, D v respectively, and their time evolution. The 
accuracy of the gray radiative transfer results can be estimated 
by comparing the local values of the deposition functions and 
net deposition to the Monte Carlo results. 

Figure 7 shows the local values of the deposition functions 
for model W7 at several times, t. The optimal value of K y at 
t ~ 60 days (t e ~ 12) for this model is K y = 0.059 Y e . This value 
of K y reproduces the deposition function in all regions of the 
model at 60 days to within a few percent. At earlier times 
(t ~ 35 days) using this value of K y there is a 6% discrepancy in 
the net deposition, D y , computed using the two methods with 
the Monte Carlo deposition function higher than the gray 
deposition function by ~6% at the peak of the distribution. At 
even earlier times, the fit is actually better due to the degener- 
acy of K y values that produce a good fit (see above). At later 
times, t = 90 days, the deposition computed using the value 
K y = 0.059^. results in a net deposition ~6% higher and local 
deposition functions less than 10% larger than those computed 
using the Monte Carlo calculations. These results are indica- 
tive of those found for other models. Figure 8 shows that 
changing K y by ~20% from the best-fit value at 60 days pro- 
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Fig. 7. — Comparison of the Monte Carlo and gray absorption calculations of the deposition function for model W7 of Nomoto et al. (1984). The dotted line in the 
top left panel shows the initial mass fraction of radioactive material, /„ d . The gray absorption calculations have all been performed with K y = 0.0597, cm 2 g“ 1 with 
7 e = 0.5, which is the optimal value of K y for t = 5 x 10 6 s after the explosion. Error bars denote the 1 a uncertainties in the Monte Carlo calculations using 200,000 
decay events (800,000 for t = 5 x 10 6 s). The long-dashed lines denote the Monte Carlo results computed assuming purely absorptive interactions. 



Fig. 8. — Deposition function for model W7 at t = 5 x 10 6 s (see Fig. 7) is 
shown for three values of the effective opacity, K y (solid lines, top to bottom : 
K y = 0.034, 0.0295, 0.025, cm 2 g _1 ). The Monte Carlo deposition function is 
also shown (dotted line). The largest percent error for k, = 0.025 or 0.034 cm 2 
g _1 is 12%. 


duces deposition functions that differ by at most 12% from the 
Monte Carlo results. 

Supernovae are initially hot and sufficiently ionized to be 
optically thick. As the supernova debris expands, the atmo- 
sphere cools and the photosphere recedes to deeper layers. 
Eventually, as the thermal diffusion timescale becomes short 
compared to the dynamic timescale, the energy deposited by 
y-rays is instantaneously balanced by thermal losses. During 
these late times an excellent approximation to the bolometric 
luminosity (not including, by convention, escaping decay 
y-rays and down-scattered X-rays) can be made by simply 
equating the instantaneous y-ray energy deposition rate to the 
bolometric luminosity. The light curve for model W7, defined 
as 

Uf) = Z k AM k d k (t)e rad ^rad D y (t)M, ai , 

is shown in Figure 9 computed from the Monte Carlo simula- 
tion and from the (time-dependent) best fit opacity gray calcu- 
lation. Light curves computed using time-independent values 
of K y = 0.02, 0.025, and 0.03 (Y e = 0.5) are also illustrated. For 
D y > 0.95, the luminosities computed using either of these four 
values for k 7 are accurate to better than ~ 3%. The light curve 
computed using the time-dependent best-fit k 7 is accurate to 
within 2% at all times. Those for K y = 0.02, 0.025, and 0.03 are 







774 


SWARTZ, SUTHERLAND, & HARKNESS 



Fig. 9. — Light curve for model W7, defined as L(r) = AM t 4 t (t)f „ d = 
e„ d Z),M„ d oc D y , is shown for the Monte Carlo simulation (solid line) and for 
the gray results using the best-fit value of k, ( dashed line ) and for three time- 
independent values of K r (0.02, 0.025, and 0.03 cm 2 g" 1 ; dotted lines). The 
luminosities computed using the best-fit value of k, differs from the Monte 
Carlo results by less than 2%. Those for K y = 0.02, 0.025, and 0.03 are in error 
by less than 22%, 10%, and 18%, respectively. The largest error occurs at 
t ~ 90 days for K y = 0.02 and 0.025 and at t = 1200 days for k, = 0.03. Only 
the “Ni and “Co decay spectrum is included. For l > 1000 days, several less 
abundant and long-lived radioactive isotopes contribute to the light curve. The 
contribution from decay positrons is not included. 

in error by less than 22%, 10%, and 18%, respectively. The 
largest error occurs at t ~ 90 days for tc y = 0.02 and 0.025 and 
at r == 1200 days for K y = 0.03. It should be noted that in the 
computation of realistic light curves one would also include the 
kinetic energy of the positrons associated with 19% of the 56 Co 
decays. If there are no radially combed magnetic field lines that 
facilitate the escape of these positrons then the positron kinetic 
energy contribution would correspond to an additional, local, 
energy deposition term ~0.04€ rad / rad which should be added 
to the results from either the Monte Carlo calculation or the 
gray absorption calculation. As noted earlier, if the magnetic 
field is radially-combed, then the positrons and most of the 
primary, high-energy electrons (produced by first and second 
scatterings of y-rays) will escape and the calculation of the 


energy deposition function will differ significantly from that 
described here. 

4. CONCLUSIONS 

The radioactive model for supernovae, in which freshly syn- 
thesized 56 Ni and its daughter nucleus, 56 Co, power the late- 
time luminosity, is well established and has been confirmed 
observationally. Numerical applications of the radioactive 
model are burdened by the complexity of the decay y-ray inter- l 

actions with the supernova material. Though physically accu- 
rate, traditional Monte Carlo techniques are excessively 
demanding of computer resources especially for light curve 
studies where the calculation of the y-ray energy deposition 
must be repeated many times. Many alternative methods have 
been invoked yet none have been rigorously justified through 
quantitative comparisons to Monte Carlo methods. 

Perhaps the simplest means of computing the y-ray energy 
deposition, which is also computationally efficient and physi- 
cally reasonable, is the gray radiative transfer method. The 
gray transfer technique was compared in detail to Monte Carlo 
simulations in this work. Solutions to the gray transfer equa- 
tions can be matched to the corresponding Monte Carlo 
results by adjusting a single parameter: the energy- 
independent material opacity, K y . 

We applied a merit function, Q, similar to a ^-statistic, to 
systematically identify the best-fit value for K y , which was 
found to be weakly dependent on the optical thickness of the 
supernova atmosphere. This dependence can be explained 
qualitatively by the fact that as photons downscatter through 
successive interactions the opacity increases with the result 
that the effective opacity is larger for greater optical depths. 

This dependence also reflects the fact that the gray radiative 
transfer model is not fully appropriate to the y-ray deposition 
problem. The accuracy of the gray transfer solution fits to the 
Monte Carlo deposition functions was also analyzed. Using 
the optimal values of K y , obtained by minimizing Q, produce 
fits to a typical Type la supernova light curve that are accurate 
to within 2% of the Monte Carlo results for times t from 
maximum light (t ~ 12 days) to the latest phase included in the 
study (t = 1200 days). Neglecting the weak optical depth 
dependence of K y by using the asymptotic value of K y = 0.05 Y e 
(the optimal value for r e < 1) also produces an acceptable fit to 
the light curve with errors not exceeding 10%. 
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